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We present a new method for isothermal rigid body simulations using the quaternion representa- 
tion and Langevin dynamics. It can be combined with the traditional Langevin or gradient (Brow- 
nian) dynamics for the translational degrees of freedom to correctly sample the NVT distribution 
in a simulation of rigid molecules. We propose simple, quasi-symplectic second-order numerical in- 
tegrators and test their performance on the TIP4P model of water. We also investigate the optimal 
choice of thermostat parameters. 



I. INTRODUCTION 

Classical molecular dynamics simulation of an isolated 
system naturally samples states from a microcanonical 
(NVE) ensemble, where the number of particles N, vol- 
ume V, and total energy of the system E are held con- 
stant. However, in many cases it is desirable to study 
the system in a more experimentally relevant canonical 
{NVT) ensemble, where the temperature T is specified 
instead of E. In order to sample from the canonical 
ensemble, the molecular dynamics equations of motion 
are modified by introducing the interaction of the sys- 
tem with a "thermostat" . There exist a large variety 
of approaches for introducing such a thermostat, which 
can be roughly classified into two categories: determin- 
istic and stochastic, depending on whether the resulting 
equations of motion contain a random component (for a 
review, see, e.g. Ref. 

Among various deterministic approaches, those based 
on coupling the system to a few external degrees of free- 
dom (e.g. Nose-Hoover thermostat) have become very 
popular. Given ergodicity in the molecular system dy- 
namics, such thermostats are proven to generate correct 
canonical ensemble sampling of the system phase space. 
However, since the thermostat variables are coupled and 
control directly only global system quantities (e.g. ki- 
netic energy), such thermostats rely on the efficient en- 
ergy transfer within the system to achieve equipartition 
within the canonical distribution, such that the average 
energy of each degree of freedom within the system is 
equal to ksT. Therefore, in a system where the energy 
transfer between its different parts is slow (e.g., systems 
combining fast and slow degrees of freedom), the simple 
Nose-Hoover thermostat may have difficulty maintain- 
ing the same temperature for the different parts of the 
system. In this case more complicated thermostats are 
necessary, for example, Nose-Hoover chain thermostat, 
or separate thermostats for different parts of the systems 
(see, e.g. Ref. |3j). 

The stochastic approach exploits ergodic stochastic dif- 
ferential equations (SDEs) with the Gibbsian (canon- 
ical ensemble) invariant measure. For this purpose, 
Langevin-type equations or gradient systems with noise 
can be used (see, e.g. Refs . l6ll8l lT2lll^ll6l and references 
therein). Stochastic thermostats, with their independent 
thermalization of each degree of freedom, provide direct 



control of equipartition and thus do not need to rely on 
the efficient energy transfer within the system. 

In order to achieve such a direct thermalization of the 
system, one needs to be able to apply stochastic ther- 
mostats to all types of degrees of freedom. The standard 
Langevin equations for translational degrees of freedom 
are well known, while Langevin thermostats for systems 
with constraints have been proposed quite recently^i?. In 
this paper we introduce Langevin equations for the rigid 
body dynamics in the quaternion representation and pro- 
pose effective second-order quasi-symplectic numerical 
integrators for their simulation. These equations can be 
coupled either with Langevin or Brownian dynamics for 
the translational degrees of freedom. 

In Section |TT] we recall the Hamiltonian system for rigid 
body dynamics in the quaternion representation from 
Ref. ,2, based on which we derive Langevin and gradient- 
Langevin thermostats. Second-order (in the weak sense) 
numerical methods for these stochastic systems are con- 
structed in Section [1111 We test the thermostats and the 
proposed numerical integrators on the TIP4P model of 
water^^. The results of our numerical experiments are 
presented in Section IIVI In particular, we investigate 
the optimal choice of thermostat parameters and the dis- 
cretization error of the numerical methods. A summary 
of the obtained results is given in Section [V] 



II. EQUATIONS OF MOTION 

We consider a system of n rigid three-dimensional 
molecules described by the center-of-mass coordinates 



„1T 



'2' '3 



and 

the rotational coordinates in the quaternion representa- 
tion q = (glT,...,g«T)T g ^4„^ ^ {qlqlqi,qlY e 

R'', such that \q^ = 1. We use standard matrix nota- 
tions, and "T" denotes transpose. Following Ref. 12, we 
write the system Hamiltonian in the form 



^^(r,P,q,7r) 



2m 



-^^^(g^^^)-^C/(r,q), (1) 

3 = 1 1=1 



where p =(p 



IT 



p3n 



= {pi,P2,P3V ^ 



in 



are the center-of-mass momenta conjugate to r, tt 



IT 



p4n 



are 



the angular momenta conjugate to q, and U (r, q) is the 
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potential interaction energy. The second term represents 
the rotational kinetic energy of the system with 

Viiq,w)^^[7r^Siq]\ g, ^ G M^ ^ = 1,2,3, (2) 

where the three constant 4-by-4 matrices Si are such that 

Siq = (-gi, go, 93, -92)"^, S2q = (-92,-93, 90, 91)"^, 
'S'39 = (-93,92,-91,90)^, 

and Ii are the principal moments of inertia of the rigid 
molecule. The Hamilton equations of motion are 

at m at 



dq- 



1=1 



(3) 



^=-E^,^^(9^^')-V,.C/(r,q), 
j = 1, . . . ,n. 

It is easy to check that if the initial conditions are chosen 
such that |9-'(0)| = 1, then the corresponding Hamilton 
equations of motion ensure that 

|9^(0I = 1, J = foraUt>0. (4) 

In the rest of this section we derive stochastic ther- 
mostats for this molecular system, which preserve (j4|. 
They take the form of ergodic SDEs with the Gibbsian 
(canonical ensemble) invariant measure possessing the 
density 

p(r, p, q, tt) cc exp(-/3i/(r, p, q, tt)), (5) 
where /3 = l/lksT) > is an inverse temperature. 

A. Langevin-type equations 



where 7 > and F > with 7r > are the fric- 
tion coefficients for the translational and rotational mo- 
tions, respectively, measured in units of inverse time, 
which control the strength of coupling of the sys- 
tem to the "heat bath" ; g is a 3-dimensional appro- 
priately normalized vector; G is a 4-dimensional vec- 
tor, which provides a balance in coupling various ro- 
tational degrees of freedom with the "heat bath"; b 
and B are 3-by-3 and 4-by-4 matrices, respectively; 

and (v^fT,wT)T ^ {w^^,...,w"'^,W^^ ,...,W"^y is 
a (3n -f 4ri)-dimensional standard Wiener process with 
= iwi,wi,wi)'^ and = {W^i, ,W( ,Wi,Wiy . 

For simplicity, we assumed here that g, G, b, and B are 
the same for all n molecules, although one could choose 
them depending on the molecule number j. The lat- 
ter can be especially useful for systems consisting of sig- 
nificantly different types of molecules. It is also natu- 
ral to require that each degree of freedom is thermal- 
ized by its own independent noise, and in what fol- 
lows we assume that the matrices b and B are diago- 
nal. Further, we suppose that the coefficients of ©-([T]) 
are sufficiently smooth functions and the process X{t) — 
{'R^{t),P'^{t),Q^{t),U^{t)y is ergodic, i.e., there exists 
a unique invariant measure fj, oi X and independently of 



p 14n 



there exists the limit 



lim E(p{X{t;x)) = / <f{x)dfi{x) := ip'"''^ (8) 



for any function ip(x) with polynomial growth at infinity 
(see Refs. IrHs HlGlflsl and references therein). Here X{t; x) 
is the solution X{t) of I©-® with the initial condition 
X{0) = X{0;x) = X. 

It is not difficult to see that the solution of ([6])- ([7]) 
preserves the property (H]), i.e., 



Consider the Langevin-type equations (in the form of 
Ito)^ 



|g^(t)|-i, j-i,. 



for aU t > 0. 



(9) 



dR^ = —dt, R^{0) = 



(6) 



dP^ = ~Vr,U(R,Q)dt 

--ig{P\W)dt + b{W)dw^{t), P^'(O) = 
3 

dQ^ = Y.^..Vi (Q^" , )dt, (0) = 9^ W\^l, 



1=1 



dff = -J2^q^ViiQ','n')dt~Vq.U{R,Q)dt (7) 



1=1 



- rG{Q^,W)dt + B{Q^,W)dW'{t), W{0) = TT^ 
j = 1, . . . ,n. 



Now we find relations between 7, F, g, G, 6, and B 
such that the invariant measure fj, is Gibbsian with the 
density ([5]). The density p(r, p,q, tt) should satisfy the 
stationary Fokker-Planck equation: 



L*p = 0, 



(10) 



where 



3 



1 9^ 
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1=1 



1=1 



After some calculations, we get the required relations: 



TO 2 



m 



+ 7 



dp>^ 



and 
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(3 , 

TO 



(11) 



dBudH_ 



(12) 



r — J - r/3G, — - = 0. 



Since numerical methods are usually simpler for sys- 
tems with additive noise, we limit computational con- 
sideration of thermostats in this paper to the case of bu 
and Bu both being constant. At the same time, we note 
that general thermostats ([6])- ([7]) with (|lip - p^ may have 
some beneficial features for certain systems but we leave 
this question for further study. For constant bu and Bu^ 
the relations PT|) - P^ take the form 



75, , ) = ^ f and TG, {q' , ] 

TO 2 



/3 



2 dnl' 



In the considered molecular model it is natural to have 
the same value for all bu, i — 1, 2, 3, and the same value 
for all Bu, i — I, . . . ,A. Further, taking into account the 
form of the Hamiltonian ([T]) and that 

3 3 

1=1 1=1 ' 

1 ^ 1 

= -^"^jSiqiSiqfn, 

1=1 ' 

we can write 

G{q,n) = J{q)n and B^ = 

with 

. ELl7;Siq[Siq]' 

J{q) = =^ 1 and M 

1^1=1 ii 

Note that Tr J(q) = jqp = 1. 



2Mr 



Y^3 1 

2^1=1 /, 



(13) 



-Vr-(pp) + Vp-[Vrt/(r,q)p] 



Thus, in this paper under 'Langevin thermostaV we 
understand the following stochastic system 



dW ^ —dt, RHQ) = r^, 

TO 

dP^ = -Vr:>U{R,Q)dt 



(14) 



2to7 



dQ' = Y.^.oVi{Q\W)dt, Q^(0) = g^ |g^| = 1,(15) 
1=1 

3 

dW - -^V5.1^(Q^ff )dt- V,. C/(R,Q)di 



1=1 



I2MT 



-TJ{Q^)Wdt+ ^^-dW^{t), ff(0) = 7r^ 
j = 1, . . . ,n, 

where J{q) and M are from (fT3)) . the rest of the notation 
is as in ([6])- ([7]). We recall that 7 and F are free parameters 
having the physical meaning of the strength of coupling 
to the heat bath. 

Let us fix a molecule and write the equations for the 
body-fixed angular velocities lu^, ujy, and tUz correspond- 
ing to the rotational Langevin subsystem To this 
end, we recall^, that 

uj, = 2iSiQ)'^Q, ujy^2{S2QVQ, Lo,^2{S:iQYQ. 

Then we obtain 

, (ti h-h \,, MF 

dijJx = 1" ^ T ^y^z ji-^xdi 

\li li I ill 



1 2A/F 



Ii\ 13 



, (t2 , /3-/1 \,. MT 
dujy = — H LJxUJz dt — —uoydt 



1 /2MF ,^ 



(16) 



dw, = — 



T3 Ii-h \ .. MT 



— H ^xLOy dt - 

h h I 



Ah 



.dt 
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where are the torques, — —^{SiQ^Wqll, and dTi = 
^ X]^=i c^W^jj which can be interpreted as random 
torques. For F = 0, (fTB)) coincide with the equations for 
the angular velocities in Ref. 2". We also note that due to 
the form of the Hamiltonian ([T]) the auxiliary velocity luq 
used in Ref. [1 in the derivation of the Hamiltonian system 
for rigid-body dynamics is identically equal to zero. 



where all the notation is as in p^ - p^ and, in particular, 
J{q) and M are from (jl3[) . The invariant measure of 
([ni)-(|IHl) is ® integrated over p. The property © is 
preserved. The gradient-Langevin thermostat has two 
free parameters, > and F > 0. The latter is the 
same as in (jl4p - (jl5p . while the former, measured in units 
of time, controls the speed of evolution of the gradient 
subsystem ([T7|. 



B. A mixture of gradient system and 
Langevin-type equation 

Another possibility of stochastic thermostating of ([3]) 
rests on a mixture of a gradient system for the transla- 
tional dynamics and Langevin-type equation for the ro- 
tational dynamics. We note that according to the density 
of Gibbsian measure ([5]) the center-of-mass momenta P 
are independent Gaussian random variables and they are 
independent of the other components of the system, so 
we can avoid simulating P via a differential equation. 

Consider the 'gradient-Langevin thermostaV 



It is important to note that the gradient system does 
not have a natural dynamical time evolution similar to 
Hamiltonian or Langevin dynamics. This is because 
changing parameter v simply leads to a time renormal- 
ization of the gradient subsystem p7p . However, when 
linked with the Langevin dynamics for rotational degrees 
of freedom, as in pT)) - p^ . parameter v controls the 
"speed" of evolution of the gradient subsystem relative 
to the speed of the rotational dynamics. 



dR 



-VrC/(R, q,)dt- 



-.dMt), R(0) = r. 



To check that 



(17) 

3 

dQ' - V,. ^ Vi{Q^,W)dt, Q^O) = q\ W\ = 1, 
1=1 

(18) 

3 

dff = -V,. Vi{Q^,W)dt - V,. C/(R, Q)dt 



n 3 



- TJ{Q3)Wdt 
j = l,...,n, 



I2MT 



dW^{t), ff(0) = 7r^ 



p(r, q, tt) a exp (-/3 ^ ,tt^) + U (r, q) 



3=1 1=1 



is the density of the invariant measure for ([T7])-([TH]), one 
needs to consider the Fokker-Planck equation (jlOp with 
the following operator: 



-V,. • {v^.Yyi{q\T:')i)j I + ^Vr • [Vrf/(r,q)p] 



1=1 



Let us remarki^ that the gradient sub-system ([TT)) can be viewed as an overdamped limit of the Langevin trans- 
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lational sub-system ^1)1 for a fixed Q. 



III. NUMERICAL INTEGRATORS 

In this section we consider effective second-order nu- 
merical metliods for tfie Langevin thermostat (|14p - (fT51) 
and the gradient-Langevin thermostat (fT7)l - (|18p . We 
first recall the idea of quasi-symplectic integrators for 
Langeviri-type equations introduced in Ref. [lO (see also 
Refs. [Tllfl2 ) and also some basic facts from stochastic 
numericsii. 

Consider the Langevin equations (fTil) - P^ . Let Dq £ 
R'', d = 14n, be a domain with finite volume. The trans- 
formation X = (r^, p^, q^, TT^)^ I— > X{t) = X{t;x) = 
CR^it; x), P~^{t; x), Cl^{t; x),U^{t; x)y maps Do into the 
domain Dt- The volume Vt of the domain Dt is equal to 



Vt = dX^ ... dX"^ 



D{X'^,...,X° 



(19) 



D{x^,...,x'i) 



dx' 



.dx" 



The Jacobian determinant J is equal to (see, e.g., Ref.[ 



DiX\...,X'') 
D{x^,...,x'i) 



exp (-n(37 + T)-t). (20) 



The system p^ - p3|) preserves phase volume when 7 = 
and r = 0. If 7 > and r > with 7r > then phase- 
volume contractivity takes place. 

If we omit the damping terms, —^P'-' and — FJII^, in 
P^ - (|15p then the system becomes a Hamiltonian sys- 
tem with additive noise^iii, i.e., its phase flow preserves 
symplectic structure. Under 7 = and F = 0, p^ - pS]) 
takes the form of the deterministic Hamiltonian system 
®. 

We say that the method based on a one-step approx- 
imation X — X{t -\- h;t,x), ft. > 0, is symplectic if X 
preserves symplectic structure^iii. It is natural to expect 
that making use of numerical methods, which are close, 
in a sense, to symplectic ones, has advantages when ap- 
plying to stochastic systems close to Hamiltonian ones. 
In Ref. [13 (see also Ref. [TT| ) numerical methods (they 
are called quasi-symplectic) for Langevin equations were 
proposed, which satisfy the two structural conditions: 

RLl. The method applied to Langevin equations degen- 
erates to a symplectic method when the Langevin 
system degenerates to a Hamiltonian one. 

RL2. The Jacobian determinant J = DX/Dx does not 
depend on x. 

The requirement RLl ensures closeness of quasi- 
symplectic integrators to the symplectic ones. As it is al- 
ways assumed, a method is convergent and, consequently. 



J is close to J at any rate. The requirement RL2 is natu- 
ral since the Jacobian J of the original system (|14p - (fT5l) 
does not depend on x. RL2 reflects the structural prop- 
erties of the system which are connected with the law of 
phase volume contractivity. It is often possible to reach 
a stronger property consisting in the equality J = J. 

We usually consider two types of numerical methods 
for SDEs: mean-square and weak^l. Mean-square meth- 
ods are useful for direct simulation of stochastic trajec- 
tories while weak methods are sufficient for evaluation of 
averages and are simpler than mean-square ones. There- 
fore, weak methods are most suitable for the purposes of 
this paper. Let us recalUi that a method X is weakly 
convergent with order p > if 



\Ecp{X{T))-Ecp{X{T))\<ChP, 



(21) 



where ft > is a time discretization step and is a suffi- 
ciently smooth function with growth at infinity not faster 
than polynomial. The constant C does not depend on ft, 
it depends on the coefficients of a simulated stochastic 
system, on ip, and T. 



A. Numerical schemes for the Langevin thermostat 

We assume that the system p^ - (IT5]) has to be solved 
on a time interval [0, T] and for simplicity we use a uni- 
form time discretization with the step ft = T/N. Using 
standard ideas of stochastic numeric o^°'^-^ including split- 
ting techniques and the numerical method from Ref. 
for the deterministic Hamiltonian system ([3]), we derive 
two quasi-symplectic integrators for the Langevin system 
dH-dUl). 

The first integrator (Langevin A) is based on split- 
ting the Langevin system p^ - p3|) into the Hamiltonian 
system with additive noise (i.e., p^ - P^ without the 
damping terms) and the deterministic system of linear 
differential equations of the form 



P = -7P 

^' - -rj(g-'")^^ j = l. 



(22) 



We construct a second-order weak quasi-symplectic in- 
tegrator for the stochastic Hamiltonian system^iiii and 
appropriately concatenat a^°i^^ it with the exact solution 
of (im. The resulting numerical method is given below. 

Introduce the mapping ^'i(t;g,7r) : (q, tt) (Q,-/7) 
defined by 



Q = cos{xit)q + sm{xit)Siq , 
n = cos(xii)7r + s\u{xit)Si-K , 



(23) 



where 



Xi 



4.L1 



Siq . 



The first quasi-symplectic scheme for (fT4] 
written in the form: 



ifTSl) can be 
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Langevin A 

Pn = p, Ro = r, Qo = q, Ho = it, 

Vi,k = Pfcexp(-7V2), 

exp(-rj(Qi)V2)ni, j = 



92 



(24) 



- 



V2± 



Rfc+1 = 



4 



Rfc + -V2,k, 
m 



(Qik^Kk) 

Mi 



(Qi.k^ni,) 
(Qik^Kk) 



^h,k 



^3{h/2;Qlni,), 
^2ih/2;Ql,,ni,), 

Mh/2;Qi,,ni,), J 



-V^3 UCRk+i, Qfc+i) 



. .,n, 



^/h 2MT , r 



3,k 



2 

V2,k - 



4 /? 



Qk+n j ~ 1' 



. ,n, 



-Vrf/(Rfe+l, Qfe+l) + —\J^^k, 



Pfe+1 = 7'3,fcCxp(-7V2), 

exp(-rj(gi+JV2)i78^;fe, J 



0, . 



where = (^i^^, . . ■,^3n,kV and 77^ = {t]{j^, , 
j = 1, . . . , n, with their components being i.i.d. with the 
same law 



P(6' = 0) = 2/3, P{9 = ±V3) = 1/6. 



(25) 



It is easy to checki^iii that the scheme is quasi- 
symplectic. Moreover, the Jacobian J of the correspond- 
ing one-step approximation is exactly equal to the Jaco- 
bian J of the original system (fT4 |) -(fT5 l) . 

To prove the second order of weak convergence of ([M]) - 
psp . we compared the corresponding one-step approxi- 
mation with the one-step approximation corresponding 
to the standard second-order weak method for SDEs with 
additive noise from Ref. MiP- 113]. The following prop- 
erties are used in this proof: 

i=o * 1=1 



dnidn 



and 



_d_ 

dx 



) = TT^(9''^')=0 forj^L 



As it is usual in stochastic numericsii, we prove conver- 
gence of a numerical method under the global Lipschitz 
assumption on the coefficients of the stochastic system, 
which can then be relaxed using the concept of rejecting 
exploding trajectoriesi^. 

Analogously to the deterministic case^, one can verify 
that the scheme ([24]) preserves ([9]), i.e., IQ-^I = 1, j = 
1, . . . , n , for all k. We summarize the properties of the 
method in the following statement. 



Proposition 1 The numerical scheme ^EJj-^E^ for 
l\14i - l\15]\ is quasi- symplectic, it preserves the structural 
property 0, and is of weak order two. 



We note that one can choose ^j. = (^ 



1 ^3n,k)~^ 



and rff^ = Wi^k^ ■ ■ ■ ^vi,kV , j = 1,...,", so that their 
components are i.i.d. Gaussian random variables with 
zero mean and unit variance. In this case the weak order 
of the scheme remains second as when we use the simple 
discrete distribution (|25p . Since simulation of the dis- 
crete random variables is cheaper than Gaussian ones, it 
is preferable to use (^5]) and it was used in all our experi- 
ments in this paper. Let us remark in passing that in the 
case of Gaussian random variables the above scheme also 
converges in the mean-square sensed with order one. 

Note that exp{—TJ{q)h/2) in is the exponent of 
a matrix. It can be computed using a standard linear 
algebra package (such as LAPACK). Since J{q) is a sym- 
metric matrix, LAPACK 's dsyev routine can be used to 
obtain the eigen decomposition 



Jiq) = Tiq)Ajiq)T^iq), 



(26) 



where T{q) is a matrix whose columns are the eigenvec- 
tors of J(g) and 

Aj{q) = diag(A,74, . . . , Aj,4) 

is a diagonal matrix of the corresponding eigenvalues. 
Then 

exp(-rj(9)/^/2) = T{q) exp(-rAj(g)V2)TT(g) , 
where 



exp(-rA,7(q)/i/2) = diag(e 



-r\j,ih/2 



-rA,7,4/i/2\ 



Alternatively, the matrix exponent exp(—TJ{q)h/2) in 
(I24p can be approximated via the Taylor expansion. To 
ensure the second-order convergence, it is sufficient to ap- 
proximate it with accuracy 0{h^) at one step; the scheme 
will remain quasi-symplectic but the Jacobian J will no 
longer be equal to JJ in ([^0]) . 
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When the parameters 7 and T are large (the strong 
couphng to the "heat bath" conditions), we propose to 
use a numerical integrator for the Langevin system JT^ 
based on the following splitting: 



dPi = -7P/ dt 



2m7 



dw(i), 



/ 2 Mr 

dU] = -TJ{q)Il]dt + J——dW\t); 



dUii = — — dt 



(27) 



dPii = - VrUiRii, Q//)di, 

3 

dQh =v,.5]^(Qi„n^„)di, 

dn^^^ = -V,.{/(R//,Q//)di 

3 

-v,.^v^(Qi„n^,,)dt, 
(=1 

j =1, . . . ,n. 

The SDEs have the exact solution: 

Pj{t) = PKO)e-^* + ^1^1^ e-^(*-^)dw(.), 

n^,(t) = exp(-rj(g)t)n^,(o) 

exp(-rj((j)(i-s))(iW^^(s). 



(28) 



(29) 



P Jo 



To construct a method based on the splitting (P7|) - (P5)) . 
we take half a step of (p7|) using ((29)l . one step of a sym- 
plectic method for ([28|) . and again half a step of ([27|) . 

The Ito integral in the expression for IIj in (|29|) is a 
four-dimensional Gaussian vector with zero mean and the 
covariance matrix 



2MT /■* 

C{t;q)^^ exp[-2rj(g)(t-s)]ds 
= ^r(g)Ac(t;g,r)rT(g), 

where T{q) is as in (^5]) and 

Ac(i;g,r) = diag(Ac,i, . • • , Ac,4) 

with 



(30) 



j 2rt , if Aj,j = 

Ac,i(^;9,r) < l-cxp(-2rA.j,.(g)t) 
I A,,,,(g) 



otherwise. (31) 



i = l,...,4. 



We note that at least one eigenvalue of J{q) equals zero 
by definition. 



Finally, introduce a 4 x 4-dimensional matrix a{t,q) 
such that 



ait;q)a'^it;q)^C{t;q). 



(32) 



Since C{t; q) is a symmetric matrix, a{t; q) can be deter- 
mined as a lower triangular matrix in the Cholesky de- 
composition of C{t;q). LAPACK's dpotrf can be used 
for this purpose. 

With the above definitions, we obtain the following 
quasi-symplectic scheme for (fT4)) - (|T5|) : 

Langevin B 
Po = p, Ro = r, Qo = q, no = TT, 



7'i,fc-Pfee-^''/2 + ^_(l_e-7'.)^ 

Kk - cxp ( - TJ{Qi)h/2)lli 

+ a{h/2;Qi)r]l, j = l,...,n, 



(33) 



V2,k =Vl,k - ^VrUiRkMk), 

nik = Kk - ^ V,. U{Rk, Qk), J = 1, . . . , n, 

R-fc+l = R-fc H 'P2,k, 



(2U.^3%) = ^3(V2;Qi,^i;,), 
iQik,nlk)-'^2{h/2;Qi,,nl,) 



(Qi+v nl,) = ^3(V2; Qi„ n^,), j = i 



P3,fc = ^2,fc-^VrC/(Rfe+i,Qfc+i), 



/5 

Tli^,=e^p{^TJiQi^,)h/2)nl, 
fc = 0, ...,7V-1, 



where ^fe (Cl,fe, • • ■ , ^3n,fe)"^. Cfc = (Cl,fe, • • ■ , C3n,fc)"^ 

and -nl {vik^y,vikV, 4 = (^i,fe>---><fe)"^> i = 
1, . . . , rt, with their components being i.i.d. with the same 
law dSS]). 

As in the case of the scheme ([M)) -(P5 |) . the Jacobian J 
of the one-step approximation corresponding to the inte- 
grator ((33| , (|25|) is exactly equal to the Jacobian J of the 
original system (fT4|) - (fT5|l . The following proposition can 
be proved. 
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Proposition 2 The numerical scheme (jggp . (jgjp /or 
([23"CI3) quasi- symplectic, it preserves the structural 
property (O, anrf is of weak order two. 

We note that if we omit the rotational component in 
([T1])-([T5]), the scheme ([11])- ([25]) coincides with a second- 
order wea k quas i-symplectic method from Ref. [13 (see 
also Refs. Illlll2[ ) and the scheme ((551) . (1^ is close to 
the one from Ref. Hi (see also Refs. 6,10Jl.il). Both 
stochastic integrators ([^ - (^5]) and (1^ degener- 

ate to the deterministic scheme from Ref. 12| when 7 = 
and r = 0. The scheme ([55)) . (^5)) is usually preferable 
when 7 and/or F are large (see our experimental results 
in Section [IV] and a discussion in the case of translational 
Langevin equations in Ref. [isl ) . It is slightly more expen- 
sive than ([M[) - (l25p due to the need of generating the ad- 
ditional 7n random variables j, and i;^ per step and 
computing Cholesky factorization. However, for most 
molecular system of practical interest in computational 
chemistry and physics, where majority of the computa- 
tional effort is spent on force calculations, the additional 
cost is negligible. 

B. Numerical sclieme for tlie gradient-Langevin 
system 

To construct the numerical scheme for the gradient- 
Langevin system ([T7 | - ([TS]) . we exploit the Runge-Kutta 
method of order two for equations with additive noise 
from Ref.lulb- 113] to simulate the "gradient" part ([T7|) 
and the "Langevin" rotational part ([T8[) is approximated 
in the same way as in ([55]). The resulting second-order 
weak scheme has the form 

gradient-Langevin 

Ro - r, Qo - q, Uo = TT, (34) 
= exp ( - rJiQi)h/2)Ui + a{h/2- Qj;)^^, 

j = 1, ... ,71, 

AR, = _|iiVrL/(R,,Q,.) + ^,&„ 
2 m 2 V "^P 

7^fe = Rfe + 2 X ARk, 
nlk = Kk - ^^9^ U{-Rk, Qfc), j = 1, . . . , n, 

(Qik^Kk) = Mh/2;Qi,nl,), 

iQik^nlk) - ^2{h/2;Q{„ni,) 
(Qi+i: nlk) = *3(V2; Qife, nl,), j = i, . . . , n, 



Rfc+i = 7^fc-Ai?, + — W— 
2 V 

-^-Vr^7(7^fc,Qfe+l), 

2 m 

Kk = ^7%-^V,.'7(Rfc+i,Q,+i), 
ni+i = exp ( - TJ{Qi^,)h/2)nl, + a{h/2- Ql+,Wk, 

j = i, fc = o,...,7V-i, 

where = (f^ fc, . . . ,<^3„,fe)"^ and r]l = {r]{,^, . . ■,^1^)'^, 

4 = ("^Ifc' ■ ■ ■ '"^Ifc)^' j = l,...,n, with their compo- 
nents being i.i.d. with the same law ([25| . The following 
proposition can be proved. 

Proposition 3 The numerical scheme <\S4l , (jgjp for 
(| j7p- (| J<$p preserves the structural property (O a?^rf is 0/ 
weafc order two. 

We draw attention to the fact that the above gradient- 
Langevin scheme requires two force calculations per step 
and thus is approximately twice as expensive as the 
Langevin schemes presented in Section IIII Al 

C. Computational errors 

Let us recall that the objective is to compute highly 
multi-dimensional integrals with respect to the Gibb- 
sian measure ^J.{x) with the density The considered 
stochastic systems ([n[) - ([T5]) and ([T7[) - (IT5[) are assumed 
to be ergodic with the Gibbsian invariant measure and 
we can represent the integrals of interest as (cf. ([8|)): 

ip'^'-a = [ ip{x)dfi{x) = lim Eip{X{t;x)). (35) 

J t^oo 

We are interested here in systems solutions of which sat- 
isfy a stronger condition, namely they are exponentially 
ergodic, i.e., for any x € R^^" and any function ip with a 
polynomial growth: 

\Eip{X{t;x))~<p''''s\<Ce-^\ t>0, (36) 

where C > and A > are some constants. In 
Refs. ^nM§ (see also references therein), one can find 
conditions under which Langevin equations are exponen- 
tially ergodic. 

It follows from ([36]) (and ([35]) ) that for any e > there 
exists To > such that for all T > Tq 

\Eip{X{T;x))-ip^^^\<e. (37) 

Then we can use the following estimate for the ergodic 
limit (^^"^f: 

(^^•^9 « Eip{X{T;x)) ^ Eif{XiT;x)) (38) 

« ^-^:=iX^^(x(')(r;x)), 
1=1 
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where T is a sufficiently large time, X is an approxima- 
tion of X, and L is the number of independent approxi- 
mate realizations. The total error 



D ^ erg ern 



(39) 



consists of three parts: the error e of the approximation 
(^£'•9 Eip{X{T] x)); the error of numerical integration 
ChP (see (HI])), and the Monte Carlo error; i.e., 

or more specifically 



0{1/L). 



Each error is controlled by its own parameter: suffi- 
ciently large T ensures smallness of the error \Lp'^^3 — 
Eif{X{T\x))\] time step h (as well as the choice of nu- 
merical method) controls the numerical integration error; 
the statistical error is regulated by choosing an appropri- 
ate number of independent trajectories L. 

The other, commonly used in molecular dynamics, nu- 
merical approach to calculating ergodic limits is based 
on the known equality 



lim — 

t— ►oo t 



ip{X{s;x))ds = if''^^ a.s., 



(40) 



where the limit does not depend on x. Then by approxi- 
mating a single trajectory, one gets the following estima- 
tor for 



erg 



ip{X{s;x))ds - (^"'■3 := 



1=1 



piX{lh;x)), 
(41) 

where T is sufHciently large and Lh — T. In Ref. |T3 this 
approach was rigorously justified in the case of ergodic 
SDEs with nondegenerate noise and globally Lipschitz 
coefficients. Let us emphasize that T in (|4ip is much 
larger than T in ((37|) and ((38|) because T should be such 
that it not just ensures the distribution of X{t) to be close 
to the invariant distribution (like it is required from T) 
but it should also guarantee smallness of variance of (f'^^^. 
See further details concerning computing ergodic limits 
in Ref. Il3 and references therein. 



of the thermostat properties on the choice of parameters 
7 and F for the Langevin system (|14p - (fT5l) and v and T for 
the gradient-Langevin system pT jl - fTSl) . as well as the de- 
pendence of the numerical discretization errors of the nu- 
merical schemes Langevin A, Langevin B, and gradient- 
Langevin on the integration step size h. As a model sys- 
tem, we use the popular TIP4P rigid model of water—. 
In order to speed up the simulations, both Lennard- Jones 
and electrostatic interactions are smoothly turned off be- 
tween 9.5 and 10 A. This truncation has minimal effect 
on the structure of liquid water, but leads to a lower 
estimated melting temperature^^ of 219 K. 

The two key requirements of a thermostat are: i) cor- 
rect sampling of phase space points distributed according 
to the Gibbs distribution at a desired thermostat tem- 
perature T, and ii) rapid relaxation of the system to the 
desired equilibrium state. The numerical accuracy of the 
sampling can be estimated by comparing the values of 
various system properties (e.g. kinetic and potential en- 
ergies, pressure) averaged over long simulation runs to 
those obtained with a much smaller step size h. 

To estimate how quickly the system relaxes to the de- 
sired equilibrium state we use the following simple ex- 
periment. A system of 2000 TIP4P water molecules is 
equilibrated at Tq = 220 K. Then the temperature of the 
thermostat is increased (instantaneously) to Ti = 270 K, 
and the run is continued until the system is equilibrated 
at the new temperature. We deliberately choose to simu- 
late the system at lower temperatures (close to the melt- 
ing temperature for this model of water), where the re- 
laxation of the system is expected to be slower. 

Assuming that the system is exponentially ergodic (see 
(|36p ). we can expect that any measured quantity A will 
relax from its equilibrium value Aq at Tq to the equilib- 
rium value Ai at Ti according to the approximate for- 
mula 

EA{t) = {A{t)) ^Ai + {Ao - Ai) exp(-t/TA) , (42) 

where ta is the characteristic relaxation time of the 
quantity A. The temperature switch occurs at t = 
and the angle brackets denote average over an ensemble 
of independent simulation runs. The subscript on ta in- 
dicates that different quantities may relax with different 
rates. The rate of system equilibration should be esti- 
mated from the maximum value of ta among the quan- 
tities of interest. 

The quantities we measure include the translational 
kinetic temperature 



rotational kinetic temperature 



(43) 



IV. NUMERICAL INVESTIGATION 

In this section we present a numerical study of the 
Langevin and gradient-Langevin thermostats derived in 
Section [Til In particular, we investigate the dependence 



j=l 1=1 

and potential energy per molecule 

n 



(44) 



(45) 
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280 




Time, ps 

FIG. 1: (Color online) Relaxation dynamics with transla- 
tional Langevin thermostat: 7 — 4.0 ps^^, F = 0. Thin 
lines show relaxation dynamics averaged over ten indepen- 
dent runs, thick lines (solid and dashed for translational and 
rotational temperatures, respectively) show the least squares 
fit to formula H42p . The estimated values of the relaxation 
times are tti, = 0.2 ps, rr„t = 1.9 ps, and tu ~ 3.6 ps. 



To illustrate the response of the system to the instan- 
taneous temperature change, we show in Fig.[T]the result 
of applying the Langevin thermostat only to the trans- 
lational degrees of freedom, i.e. F = in (fT4)l - (fT5l) . As 
expected, the translational kinetic temperature quickly 
relaxes to the new temperature, while the rotational ki- 
netic temperature and potential energy lag behind. To 
estimate the relaxation rates of the measured quantities, 
we use the least squares fit of the exponential function 
to the average measured quantity {A{t)). 

In all the simulations we performed, the potential en- 
ergy relaxation time is larger than that for either of the 
kinetic temperatures. Therefore, we determine the re- 
laxation time of the system to the new equilibrium state 
based on the value of tu. 

Varying the value of the translational Langevin pa- 
rameter 7, we observe that relaxation is slower for both 
small and large values of 7, with the fastest relaxation 
around 7 = 4.0 ps~^. The existence of an optimal value 
for the choice of the thermostat parameter is consistent 
with observations in Ref. [H and can be understood in 
terms of the interaction of the system with the thermo- 
stat. For small values of 7, the relaxation of the system 
is slow due to the limited heat flux between the system 
and the thermostat. For large 7, even though the ki- 
netic temperature relaxes very quickly, the relaxation of 
the configurational state of the system is apparently hin- 
dered by the disruptive influence of the random force on 
the Hamiltonian dynamics which is driving the system to 
the new equilibrium. 




10° 10' 10^ 

7, ps-i 



FIG. 2: (Color online) Translational temperature relaxation 
time for the Langevin thermostat (|14p - (ll5l ) as a function of 
the thermostat parameters 7 and F. 




FIG. 3: (Color online) Rotational temperature relaxation 
time for the Langevin thermostat (|14|) - (|15| ) as a function of 
the thermostat parameters 7 and F. 

A. Langevin Thermostats 

Next, we investigate the dependence of the system re- 
laxation time Tu on both 7 and F in P^ - fTS]) . In order 
to minimize the influence of the numerical discretization 
error, we use a relatively small time step of 0.2 fs. With 
such a small time step, the difference between Langevin 
A and Langevin B is negligible compared to the sam- 
pling error. To produce the results reported below, we 
use Langevin A. We evaluate tu on a logarithmic grid 
of 7 and F values using five independent runs at each 
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FIG. 4: (Color online) Potential energy relaxation time for the 
Langevin thermostat (|14p - p5p as a function of the thermostat 
parameters 7 and F. 
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FIG. 5: (Color online) Relaxation dynamics with 'optimal' 
choice of Langevin thermostat parameters: 7 = 4.0 ps~^, 
F = 10.0 ps~^. Thin lines show relaxation dynamics averaged 
over ten independent runs, thick lines (solid and dashed for 
translational and rotational temperatures, respectively) show 
the least squares fit to formula (|42|) . The estimated values of 
the relaxation times are TTt, = 0.28 ps, t%.^^ = 0.26 ps, and 
Tu = 2.0 ps. 



point. The results are shown in Figs. [H [31 and [D As 
expected, the relaxation speed of translational and rota- 
tional temperatures uniformly increases with increasing 
values of 7 and F, respectively. At the same time, the 
relaxation speed of the potential energy exhibits nonuni- 
form dependence on the thermostat parameters. As can 



FIG. 6: Dependence of the approximated average properties 
of a system of 2000 TIP4P water molecules on the integration 
time step h for Langevin A and B. The system is equilibrated 
with the thermostat parameters 7 = 4.0 ps~^, F = 10.0 ps~^, 
and T = 270 K. The quantities {Ttr)h, {%ot)h, and {U)h are 
denoted by circles, triangles, and squares, respectively. Solid 
and open symbols refer to Langevin A and B, respectively. 



be seen in Fig. [H the fastest relaxation of the system 
is achieved when the Langevin thermostat is applied 
to both translational and rotational degrees of freedom, 
with 7 = 2 — 8 ps^^ and F = 3 — 40 ps""'^. The relaxation 
dynamics of the system with 'optimal' choice of Langevin 
thermostat parameters is demonstrated in Fig.O In this 
case TTt, = 0.28 ps, rr,„t = 0.26 ps, and tu = 2.0 ps, 
which shows that the system relaxation is almost twice 
as fast as when the Langevin thermostat is applied only 
to translational degrees of freedom. Note that the re- 
sults shown in Fig. |4] for 7 larger than about 100 ps^^ 
or F larger than about 1000 ps~^ are not reliable due 
to excessive coupling of the system to the thermostat, 
which disrupts the Hamiltonian flow of the system. In 
this case, the relaxation dynamics is poorly represented 
by the exponential function ([42]) and thus the fits produce 
misleading values for the system relaxation time. 

Now we look at performance of the numerical integra- 
tors proposed in Section IIII Al for the Langevin thermo- 
stat ([T? | - p^ . Since both Langevin A and B are second- 
order methods, the calculated average quantities for sim- 
ulations with step size h should have the form^^'^^ 



{A)h = {A)o + CAh^ + 0{h^), 



(46) 



where (A) ^ denotes the average value of dynamical quan- 
tity A{t) calculated over a numerical trajectory with time 
step h. The dependence of {Tt^)h, {%ot)h, and (JJ)h on 
h for both Langevin A and B is illustrated in Fig. [6l It 
appears that Langevin A has larger discretization error 
for the rotational temperature and smaller error for the 
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TABLE I: Values of the coefficients Ca in the discretization 
errors (|46|) for the measured quantities with Langevin A and 
B. The system's parameters are as in Fig. [B] 





Langevin A 


Langevin B 


Crt, , K/fs^ 


-0.13 


0.06 


CT,,t , K/fs^ 


-0.85 


-0.14 


Cu , kcal/mol/fs^ 


-0.0007 


0.0056 



Langevin A Langevin B 




7, ps 

FIG. 7; (Color online) Dependence of {7tr)o and Crtr on j 
and r for Langevin A and B thermostats. 

potential energy than Langevin B. The linear dependence 
of the measured quantities on is maintained up to a 
relatively large time step of about h — 7 is. The values 
of the slopes Ca are listed in Table [D Both methods 
become unstable at about h — lOfs. 

We have investigated the dependence of {A)q and Ca 
on the thermostat parameters 7 and T, by running sim- 
ulations with time steps h — 2is and 3 fs and estimating 
these quantities from the straight line fit with respect 
to h^. We consider the fit to be justified if the higher 
order terms in (|46p are small, i.e., when the quantities 
(7^r)o and (7^ot)o determined from are equal to the 
thermostat temperature parameter T — 270 K. In Fig. [7| 
we show results for the translational temperature mea- 
surements in simulations with Langevin A and B. As ex- 
pected, {Ttr)o converges to 270 K. We note in passing 
that smaller statistical errors are observed at larger val- 
ues of 7. The behavior of Crt, for Langevin A exhibits 
a plateau for small and moderate values of 7 and T and 
then changes rapidly at values that are 'too large' for 
this system. By contrast, the translational temperature 
discretization error of Langevin B thermostat exhibits 
consistent behavior for all values of 7 and T. Similar 
differences between Langevin A and B can be seen in 
the measurements of rotational temperature and poten- 



Langevin A Langevin B 




7, ps 

FIG. 8: (Color online) Dependence of {7^ot)o and CT^at on 7 
and r for Langevin A and B thermostats. 



Langevin A Langevin B 




FIG. 9: (Color online) Dependence of (W)o and Cu on 7 and 
r for Langevin A and B thermostats. 

tial energy shown in Figs. [8] and [9l respectively. As can 
be seen from the plot of (7^ot)o in Fig. [51 the straight line 
fit also breaks down at large F values in Langevin A. 



B. Gradient-Langevin thermostat 

Here we describe numerical experiments with the 
gradient-Langevin thermostat P7| -(|18 p introduced in 
Section IIIBI Since the gradient system for the trans- 
lational motion does not include linear momenta P, the 
translational kinetic temperature 7tr is not available for 
measurement in this case. Therefore, in our numerical 
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FIG. 10: (Color online) Rotational temperature relaxation 
time for the gradient-Langevin thermostat (I17P - (|18( ) as a func- 
tion of the thermostat parameters v and T. 
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FIG. 11: (Color online) Potential energy relaxation time for 
the gradient-Langevin thermostat (|17p - (|18|) as a function of 
the thermostat parameters v and F. 



experiments we measure the rotational temperature %ot 
and potential energy per particle U as defined by ((44l) 
and ((15|) . respectively. 

For the particular system studied here, the gradient- 
Langevin numerical scheme (p4l) . (l25l) from Section [III Bl 
becomes unstable when the product hv is larger than 
about 200 fs^. Therefore, with the step size oi h = 0.2 fs 
used in our simulations, we can study the properties of 
the gradient-Langevin scheme with v up to about 1000 fs. 

We conducted the relaxation experiment, where we 
monitored {%ot{t)) and {U{t)) while the thermostat tem- 




FIG. 12: (Color online) Relaxation dynamics with gradient- 
Langevin thermostat: = 100 fs, F = 0. Thin lines show 
relaxation dynamics averaged over seven independent runs, 
thick lines show the least squares fit to formula (|42p . 



perature parameter was switched from 220 K to 270 K. 
The relaxation times tt;^^ and tu for these quantities 
were calculated for different values of v and F. The re- 
sults are shown in Figs. [TUlandfTT] As expected, the re- 
laxation time for rotational temperature decreases with 
increasing value of F. The somewhat surprising finding 
of this experiment is that, even for small values of F the 
relaxation time is much smaller here than in the case 
of Langevin thermostat (see Fig. [3]). As an illustration, 
we show the relaxation experiment for = 100 fs and 
F = in Fig. [121 The estimated relaxation time for ro- 
tational temperature, Tj-^t = 0.14 ps, is much smaller 
then the corresponding quantity for the Langevin ther- 
mostat, even though the relaxation time for the potential 
energy, = 2.7 ps, is similar. This indicates a very effi- 
cient heat transfer between the gradient dynamics of the 
translational motion and the rotational motion. 

The dependence of tu on the thermostat parameters 
for the gradient-Langevin system shown in Fig. [TT] is 
markedly different than for the Langevin system. In par- 
ticular, the relaxation time decreases with increasing v 
without reaching a minimum value within the range of 
v values explored. At the same time, there is little de- 
pendence on F, except for very large F where, similar 
to the Langevin system, the measurements of the relax- 
ation time are not reliable. Also, note that the relax- 
ation times are much smaller, reaching as low as 0.7 ps 
for V = 1000 fs, compared to the minimum value of about 
2.0 ps for the Langevin system. 

Of course, the direct comparison between the relax- 
ation speeds of gradient-Langevin and Langevin dynam- 
ics has to be taken with caution, since, as we mentioned 
in Section III Bl the gradient dynamics does not have a 
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FIG. 13: Dependence of the approximated average properties 
of a system of 2000 TIP4P water molecules on the integration 
time step h for the gradient-Langevin numerical method. The 
system is equilibrated with the thermostat parameters v = 
200fs, r = 5.0ps~\ and T = 270K. The quantities {Trot)h, 
and {U)h are denoted by circles and squares, respectively. 
Error bars reflect 95% confldence intervals in the obtained 
results estimated from block averages. 



natural evolution time. In particular, the mass of the 
molecule m does not have a specific meaning in the gra- 
dient system, since it can be rescaled to any value to- 
gether with h and v (see ((T7)) V Computationally, the 
relaxation speed depends on the time step h and, while 
with v — 1000 fs the gradient Langevin scheme becomes 
unstable for h larger than 0.2 fs, the Langevin scheme re- 
mains stable up to about h — 10 fs for the optimal values 
of 7 = 4.0ps-i and T = 10 ps^^. 

The dependence of discretization error in measured 
quantities on h for the gradient-Langevin scheme is 
shown in Fig. [T31 As in the case of Langevin A and B, we 
clearly see the linear dependence of {%ot)h, and (JA)h on 
/i^. The estimated slopes in (gH]) are Cr„t = -0.38 K/fs^ 
and Cu = —0.029 kcal/mol/fs^. Unfortunately, for this 
value of V the gradient-Langevin numerical integrator be- 
comes unstable for /i > 1 fs, which is rather small, given 
that the Langevin A and B integrator are stable for h up 
to about 10 fs. Still, given the observed efficient heat 
transfer from the gradient subsystem for translational 
dynamics to the rotational dynamics (see Fig. [T^] and re- 
lated discussion) , it might be of interest to construct nu- 
merical methods for the gradient-Langevin system (jl7p - 
PS)) with better stability properties than those of ([M)) . 



(PS)) : this has not been considered in this paper. 

V. SUMMARY 

The new stochastic thermostats presented in this pa- 
per are appropriate for quaternion-based rigid body mod- 
els. They are written in the form of Langevin equations 
and gradient-Langevin system (gradient subsystem for 
the translational degrees of freedom and Langevin sub- 
system for the rotational degrees of freedom). The ob- 
tained stochastic systems preserve the unit length of the 
rotational coordinates in the quaternion representation 
of the rigid-body dynamics. The thermostats allow to 
couple both translational and rotational degrees of free- 
dom to the "heat bath" . As it is shown in the numerical 
tests with the TIP4P rigid model of water, the Langevin 
thermostat relaxes to an equilibrium faster when not only 
translational degrees of freedom but also rotational ones 
are thermostated. It turns out that there is an optimal 
range of the strength of couphng to the "heat bath" . In 
contrast, the gradient-Langevin thermostat has a mono- 
tone dependence of relaxation time on the thermostat 
parameters. In the case of the Langevin thermostat, two 
quasi-symplectic second-order (in the weak sense) inte- 
grators are constructed and compared in the numerical 
tests. For the gradient-Langevin thermostat, a Runge- 
Kutta second-order method is proposed. All the meth- 
ods preserve the unit length of the rotational coordinates. 
The numerical experiments demonstrate the efficiency of 
the proposed thermostating technique. 

Relaxation times for the gradient-Langevin thermostat 
are smaller than for the Langevin thermostat. However, 
the numerical methods proposed for the Langevin sys- 
tem have better stability properties than the scheme used 
for numerical integration of the gradient-Langevin sys- 
tem. In our experimental study, the use of the Langevin 
thermostat together with the quasi-symplectic integra- 
tors was computationally significantly more efficient than 
thermostating via the gradient-Langevin system and the 
numerical scheme for it. 
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